# Required packages ####
#install.packages("quantreg")
#install.packages("doParallel")
#install.packages("snow")
require(quantreg)
require(doParallel)

# Folders ####
directory<-Sys.getenv("US_Ineq_Repl")
MMwd<-file.path(directory, "Scripts/MM")
Datawd<-file.path(directory, "Processed")

# Data ####
setwd(Datawd)
load("QR_mod.RData")

# Functions #### 
setwd(MMwd)
source("MMfun.R")

# Code begins ####
meth<-"pfn"

no_cores <- round(0.8*detectCores()) 
#no_cores = 64

m<-4500

# Parallel computing of QR Coef. for a set of uniform random variables [0,1]  ####
## t=0 ####
start_time <- Sys.time()
CoefRq0<-ParCoefRq(m=m,formula=formula0,d=D$Z0,w="orgwgt",mth=meth,seed=99,nodes=no_cores,ParRq=ParRq)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
# time = 54.33 mins

## t=1 ####
start_time <- Sys.time()
CoefRq1<-ParCoefRq(m=m,formula=formula1,d=D$Z1,w="orgwgt",mth=meth,seed=99,nodes=no_cores,ParRq=ParRq)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
# time = 54.42 mins

# Estimates of Wages ####
#Some W0e may be NA because some Z's may be NA
W0e<-SampW(D$Z0,CoefRq0,w="orgwgt",seed=15)$We
W0<-D$Z0[,"lrwage3"]
  

#Some W1e may be NA because some Z's may be NA
W1e<-SampW(D$Z1,CoefRq1,w="orgwgt",seed=16)$We
W1<-D$Z1[,"lrwage3"]

## Counterfactual ####
#Some W10e may be NA because some Z's may be NA
W10e<-SampW(D$Z0,CoefRq1,w="orgwgt",seed=17)$We

(max(W0,W0e,W1,W1e,W10e,na.rm = TRUE))

###
# Changes in Covariates ####
###

d0e<-density(W0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)
d0<-density(W0$lrwage3,bw="nrd",kernel="gaussian",na.rm=TRUE,
            weights=D$Z0$orgwgt/sum(D$Z0$orgwgt),from=0,to=6)
d1e<-density(W1e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)
d1<-density(W1$lrwage3,bw="nrd",kernel="gaussian",na.rm=TRUE,
            weights=D$Z1$orgwgt/sum(D$Z1$orgwgt),from=0,to=6)
d10e<-density(W10e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)


## Simulated Wages ####
SW1e<-SampW(D$Z1,CoefRq1,w="orgwgt",seed=18)
m<-4500
### Contribution of individual covariate ####
### Impact of Unionization ####
#I1=index of union and nonunion in the random sample implied by the model for t=1
I1<-list(i1=SW1e$Ze["union"]==1,i2=SW1e$Ze["union"]!=1)
#I0=index of union and nonunion in the data when t=0
I0<-list(i1=D$Z0[,"union"]==1,i2=D$Z0[,"union"]!=1)
seeds=100+1:length(I1)

W1union0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1union0e<-density(W1union0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Public Sector ####
I1<-list(i1=SW1e$Ze["pubsect"]==1,i2=SW1e$Ze["pubsect"]!=1)
I0<-list(i1=D$Z0[,"pubsect"]==1,i2=D$Z0[,"pubsect"]!=1)
seeds=200+1:length(I1)

W1pubsec0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1pubsec0e<-density(W1union0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Manufacturing Sector ####
I1<-list(i1=SW1e$Ze["manuf"]==1,i2=SW1e$Ze["manuf"]!=1)
I0<-list(i1=D$Z0[,"manuf"]==1,i2=D$Z0[,"manuf"]!=1)
seeds=300+1:length(I1)

W1manuf0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1manuf0e<-density(W1union0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Race ####
I1<-list(i1=SW1e$Ze["nonwhite"]==1,i2=SW1e$Ze["nonwhite"]!=1)
I0<-list(i1=D$Z0[,"nonwhite"]==1,i2=D$Z0[,"nonwhite"]!=1)
seeds=400+1:length(I1)

W1race0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1race0e<-density(W1race0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Gender ####
I1<-list(i1=SW1e$Ze["female"]==0,i2=SW1e$Ze["female"]!=0)
I0<-list(i1=D$Z0[,"female"]==0,i2=D$Z0[,"female"]!=0)
seeds=500+1:length(I1)

W1sex0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1sex0e<-density(W1sex0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#### Impact of Part Time ####
I1<-list(i1=SW1e$Ze["partte"]==0,i2=SW1e$Ze["partte"]!=0)
I0<-list(i1=D$Z0[,"partte"]==0,i2=D$Z0[,"partte"]!=0)
seeds=600+1:length(I1)

W1part0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1part0e<-density(W1sex0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Married ####
I1<-list(i1=SW1e$Ze["married"]==0,i2=SW1e$Ze["married"]!=0)
I0<-list(i1=D$Z0[,"married"]==0,i2=D$Z0[,"married"]!=0)
seeds=700+1:length(I1)

W1married0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1married0e<-density(W1sex0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Urban ####
I1<-list(i1=SW1e$Ze["smsa"]==0,i2=SW1e$Ze["smsa"]!=0)
I0<-list(i1=D$Z0[,"smsa"]==0,i2=D$Z0[,"smsa"]!=0)
seeds=800+1:length(I1)

W1urban0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1urban0e<-density(W1sex0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Experience ####
I1<-list(i1=(0<=SW1e$Ze["exper"])&(SW1e$Ze["exper"]<10),i2=(10<=SW1e$Ze["exper"])&(SW1e$Ze["exper"]<20),i3=(20<=SW1e$Ze["exper"])&(SW1e$Ze["exper"]<30),i4=(30<=SW1e$Ze["exper"]))
I0<-list(i1=(0<=D$Z0["exper"])&(D$Z0["exper"]<10),i2=(10<=D$Z0["exper"])&(D$Z0["exper"]<20),i3=(20<=D$Z0["exper"])&(D$Z0["exper"]<30),i4=(30<=D$Z0["exper"]))
seeds=900+1:length(I1)

W1exper0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1exper0e<-density(W1exper0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

### Impact of Education ####
#### High School ####
I1<-list(i1=SW1e$Ze[,"educ2"]==1,i2=SW1e$Ze[,"educ2"]!=1)
I0<-list(i1=D$Z0[,"educ2"]==1,i2=D$Z0[,"educ2"]!=1)
seeds=1000+1:length(I1)

W1hsch0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1hsch0e<-density(W1hsch0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#### Some College ####
I1<-list(i1=SW1e$Ze[,"educ3"]==1,i2=SW1e$Ze[,"educ3"]!=1)
I0<-list(i1=D$Z0[,"educ3"]==1,i2=D$Z0[,"educ3"]!=1)
seeds=1100+1:length(I1)

W1sobc0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1sobc0e<-density(W1sobc0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#### Associate Degree ####
I1<-list(i1=SW1e$Ze[,"educ4"]==1,i2=SW1e$Ze[,"educ4"]!=1)
I0<-list(i1=D$Z0[,"educ4"]==1,i2=D$Z0[,"educ4"]!=1)
seeds=1200+1:length(I1)

W1ad0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1ad0e<-density(W1sobc0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#### Bachelor ####
I1<-list(i1=SW1e$Ze[,"educ5"]==1,i2=SW1e$Ze[,"educ5"]!=1)
I0<-list(i1=D$Z0[,"educ5"]==1,i2=D$Z0[,"educ5"]!=1)
seeds=1300+1:length(I1)

W1bche0e<-SampWcf(SW1e,I1,I0,seeds,m)
d1bche0e<-density(W1bche0e,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

###
### Changes in Coefficients ####
###

#Unionization
We1union0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"union"),w="orgwgt",seed=18)$We
de1union0<-density(We1union0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Public Sector
We1pubsec0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"pubsect"),w="orgwgt",seed=19)$We
de1pubsec0<-density(We1pubsec0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Manufacturing Sector
We1manuf0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"manuf"),w="orgwgt",seed=20)$We
de1manuf0<-density(We1manuf0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Race
We1race0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"nonwhite"),w="orgwgt",seed=21)$We
de1race0<-density(We1race0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Gender
We1sex0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"female"),w="orgwgt",seed=22)$We
de1sex0<-density(We1sex0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Part Time
We1part0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"partte"),w="orgwgt",seed=23)$We
de1part0<-density(We1part0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Married
We1married0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"married"),w="orgwgt",seed=24)$We
de1married0<-density(We1married0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Urban
We1urban0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"smsa"),w="orgwgt",seed=25)$We
de1urban0<-density(We1urban0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Experience
We1exper0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"exper"),w="orgwgt",seed=26)$We
de1exper0<-density(We1exper0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Education
##High School
We1high0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"educ2"),w="orgwgt",seed=27)$We
de1high0<-density(We1high0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Some Bachelor
We1soco0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"educ3"),w="orgwgt",seed=28)$We
de1soco0<-density(We1soco0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#Associate Degree
We1ad0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"educ4"),w="orgwgt",seed=29)$We
de1ad0<-density(We1ad0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)

#College
We1coll0<-SampW(D$Z0,Ind_coef(CoefRq1,CoefRq0,"educ5"),w="orgwgt",seed=30)$We
de1coll0<-density(We1coll0,bw="bcv",kernel="gaussian",na.rm=TRUE,from=0,to=6)


rm(D,model0,smodel0,model1,smodel1)

setwd(Datawd)
save.image("MM_fig.RData")

